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Abstract This paper reports on a simple pure numerical method developed for computing 
Hansen coefficients by using recursive harmonic analysis technique. The precision criteria 
of the computations are very satisfactory and provide materials for computing Hansen's 
and Hansen's like expansions, also to check the accuracy of some existing algorithms. 
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1 INTRODUCTION 

Hansen coefficient Cefola 1 19771 1 is an important class of functions which frequently occur in many 
branches of Celestial Mechanics such as planetary theory(Newcomb [l 8951 > and artificial satellite mo- 
tion (Allar fT967l Hughe jT977| ). Moreover , there are extensive forms of Hansen like expansions (Klioner 
et. al. 119981 : Sharaf[l985 1986) which play important roles in the expansion theories of elliptic motion. 
Giacalia (1976) noted that Hansen's coefficients appears in satellite theory in expression of the disturb- 
ing function due to the primary and due to the presence of a third body and they are usually called 
Eccentricity Functions. He derived recurrence relation for these functions and their derivatives, as they 
appear in the evaluation of geopotential and third body perturbations of an artificial satellite. Also in 
119871 he proved Hansen's coefficients for Fourier series in terms of the mean anomaly correspond to 
a rotation of the orbital plane proportional to the eccentricity of the orbit. They are given in terms of 
Bessel functions and generalized associated Legendre functions which arise through the transformation 
of spherical harmonics under rotation. In |1981| Hughes computed tables of analytical expressions for 
the Hansen coefficients X™' ±m (e) and X~ (n+1) ' ±m (e) when 1 < n < 30 and < m < n. In [19901 
Branham derived a recursive calculation of Hansen coefficients which are used in expansions of elliptic 
motion by three methods: Tisserand's method, Von Zeipel-Andoyer method with explicit representation 
of the polynomials required to compute the Hansen coefficients and von Zeipel-Andoyer method with 
the value of the polynomials calculated recursively. Vakhidov (2000!) studied in detail efficient approx- 
imations of Hansen coefficients using polynomials in terms of the eccentricity. He and Zhang d 19901 ) 
used Hansen coefficients to compute general perturbations of the asteroids of Flora group due to Jupiter. 
Breiter et.al ( 120041 1 show that most of the theory of Hansen coefficients remains valid for X^, when 7 
is a real number, also, the generalized coefficients can be applied in a variety of perturbed problems that 
involve some drag effects. Sadov (2008) deals analytically with the properties of Hansen's coefficients 
in the theory of elliptic motion considered as functions of the parameter 77 = \/l — e 2 where e is the 
eccentricity. 

In the present paper, we develop a simple pure numerical method for computing Hansen coefficients by 



* E-mail: hassanselim@hotmail . com 



2 



M. A.Sharaf&H. H. Selim 



using recursive harmonic analysis technique. The precision criteria of the computations are very satis- 
factory. The importance of the method is that it not only provides materials for computing Hansen's and 
also Hansen's like expansions but also, it can be used due to its simplicity and accuracy, to check the 
accuracy of the different algorithms already existing. 



2 BASIC FORMULATIONS 
2.1 Properties of Least -Squares 

Let y be represented by the general linear expression of the form Ylf=i c i4>( x ) where cfi's are linear 
independent functions of x. Let c be the vector of the exact values of the c's coefficients and c be 
the least -squares estimators of c obtained from the solution of the normal equations Gc = b . The 
coefficient matrix G(L x L) is symmetric positive definite, that is, all its eigenvalues Vi; i = 1,2, .., L 
are positive. Let E C(z) denote the expectation of z and a 2 the variance of the fit, defined as 



a 1 = q n /(N - L) (1) 



where 



q n = (y-$ J c) 7 (y-$ J c) (2) 

N is the number of observations, y is a vector with elements y^ and $(L x N) has elements 
<fiik = f i{xk) ■ The transpose of a vector or a matrix is indicated by the superscript 'T'. 
According to the least- squares criterion, it could be shown that(Sharaf et.al.2000) 

- The estimators c given by the least- squares method give the minimum of q„ . 

- The estimators c of the coefficients c, obtained by the method of least-squares, are unbiased; 
i.e.EC(c) = c 

- The variance-covariance matrix Var(c) of the unbiased estimators c is given by 



where G 1 is the inverse of G. 



Var{c) = cr 2 G-\ (3) 
- The average squared distance between c and c is 

EC(D*)=a*J2± (4) 
i Vl 

2.2 Harmonic Analysis of a Periodic Function 

Let it be required to find a sum 



E 



cij cos jx + bj sinjx (5) 

which furnishes the best possible representation of a function u{x), when we are given that u(x) 
takes the values u 0) ui, iti_i when x takes x , Xi, respectively , m being some number 
greater than 2s. The problem is to determine the (2s + 1) constants, a , aj and bj;j = 1, 2, s so as 
to make the expression (5) takes, as nearly as possible, the I values u ,ui, when x takes the 

values x ,xi, Xi-\. To do so we shall make use of the method of least squares and we get 
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TjaoVoi + ajfjij + ^2 bjfiij = di, i = 0, 1, s; 

j'=i i=i 
! f s 



i=i j'=i 



where 



i-1 

Vij = Vji = /] cos ixk, i = 0, 1, •■■ , 8,j = 0, 1, —s; 

fc=0 

i-1 



/3 9J - = y^cosjxfcsingxfc.j = 0,l,-,s,g = l,2,...,s; 

i-1 

7gj = Ijq = 'Y^smqx k smjx k ,q = 1,2, ...,s,j ; = 1,2,..., a; 

fc=0 

i-1 

di = ^ Mfc cosix kl i = 0, 1, s; 

k=0 
i-1 

c i = ^2 u kSinqx kl q = l,2,...,s. (7) 

Equations (7) are the normal equations. These equations represent a set of linear equations in (2s + 
1) unknowns a's and b's coefficients and could be solved by any of the methods adopted for linear 
systems. However, the coefficient matrix of this set could be reduced to a diagonal one by certain choice 
of the arguments x k and in this case the a's and b's are determined exactly and the problem is known as 
harmonic analysis. 

In the method of harmonic analysis, the arguments x k take the special values; 

_.(,-.,.*. 

For these values, the T]'s, (3's and j's of Equations (7) become: 
For i = j ± : ^ = j jt = \l; = 0. 
For i ^ j : ifrj = 7y = Pij = 

Consequently the a's and b's coefficients could then be computed exactly from 



u k cos j. — k,j = 0,l,..,s; 



2 ^ 2?r 

b q = - 2_^u k smq. — k,q = 1,2, ...,s. (9) 

fc=0 



where = 1 if j = ; fj, = 2 if j > 
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2.3 Hansen Coefficients 

Consider elliptic motion expansions of (r/a) n cos mv and (r/a) n sin mv in terms of the mean anomaly 
M that is, 

J2Al' m cos kM; 

k=0 

^5™' m sinfcM. (10) 

fc=i 

where a is the semi major -axis , r the radial distance, n is a positive or negative integer , while m is 
positive integer and v the true anomaly in elliptic motion . The A's and B's coefficients called Hansen's 
coefficients ,are functions of the eccentricity e. 

The relations between the eccentric anomaly E and the anomalies M, v are given for elliptic motion as 
follows: 

- The relation between E and M is well know Kepler's equation 




M = E-e&mE. (11) 
- The fundamental relations between v and E in an elliptic orbit are 



v 1 + e E 
tan 2 = V T~~e tan ~2 

These equations are the most useful relations between v(E) and E(v) , since ~ and are always 
in the same quadrant. There is a possibility of numerical trouble when Equation (12) is used with 
angles that are near ± f as the two tangents become infinite. In order to avoid this difficulty, Broucke 
and Cefola [1973 1 established the formula 

1/ t-a BsvclE „ on 

tan -(«-£) = — , (13) 

2 V ' l-pcosE' V 

where 



e 1 + v/1 - e 2 

Equation (13) is free of numerical trouble, no matter what the values of the angles are. Moreover, it 
can be easily used because the angle (v — E)/2 is always less than -| for all elliptic orbits. 

- Finally the relation between r and E is 



(-) = l-ecos£ (15) 
a 

3 COMPUTATIONAL DEVELOPMENTS 

3.1 Practical Computations of the a's and b's Coefficients 

The a's and b's coefficients of Equations (9) could be computed efficiently (Ralston & 
Rabinowit j 1 9781 lfrom 

ctj = j{u + F hj cos^- j - F 2 ,j};j = 0, l,...,s, (16) 
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I ' 1 

where, for any j the F's are computed recursively from 



b q = -F hq sill — q; q = 1, 2, s. (17) 



i*k,j = u k + 2 cos XjFk+ij - F k+2 ,j (18) 
by using the initial conditions Fjj = i'i+ij = 0, starting with fc = i — 1 to compute successively 

Fi—l^j , Fi-2,j ; •• •) -^1 J ■ 

3.2 Error Estimates 

- The variance of the fit (Equation (1)) is given by 

2 S 2 

= T^l (19) 

where the sum of the squares of the residuals 6 2 is given as (Ralston & Rabinowit j 19781 1 

i-l 

I, 2 - 

2 1 

i=0 j = l 

Clearly both a and S depend on the number s of the a's and b's coefficients. If the precision is 
measured by probable error PE, then 

PE = 0.6745ct (21) 

- Since the coefficient matrix G of the harmonic analysis is diagonal with elements of the same value 
1/2 , then according to Equation (3) the standard error of each of the a's and b's coefficients is 



5 2 = Y,n 2 -\[2a 2 + j^{a 2 + b 2 )] (20) 



a coeff = a V 1 ^2) 

The corresponding probable error for each coefficient is 

PE coeff = 0.6745a coeff (23) 

- The average squared distance between the exact and the least -squares values (Equation (9)) is given 
according to Equation (4) as 

Q = CE(D 2 ) = j<j 2 (24) 

3.3 Choosing the Number of the Coefficients 

In practice , since we do not know s , we would evaluate a's and b's coefficients for s = 1,2, , then 
compute S 2 (Equation(19)), and continue as long as S 2 decreases significantly (within a given tolerance 
Tol) with increasing s . 

3.4 The Special Values 

The special values of the left hand sides of Equation (10) are computed as follows: 
1. Mi = * = 0, 1, I - 1 . 
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2. For each M solve Kepler's equation (Equation(l 1)) by Newton-Raphson iterative method (or any 
other method). Let E a be an initial approximation of E; define for k = 0, 1, 2, ... 



Each Ek + i should approximate E more closely than E^ . For the initial approximation E a use the 
value (Battii fl999l > 



e 2 < El and | H(E i+1) |< 100e 1 ,e 2 =| | if | El+l |> l;e 2 =| E l+1 - E, | if 



| Ei+i |< 1, where e\ is a given tolerance and H(E) = M — E — e sin E. 

3. For each E compute v using Equation (14) and (-)" from Equation (10) 

4. For each v compute cos(mw). 

5. Finally, find the product of the values of (^)™(of step 3)and cos(mu) (of step 4). 
3.5 Numerical Results 

The above computational developments are applied for calculating Hansen's coefficients of Equation 
(10) with input constants as I — 100, Tol = 10~ 6 and E\ — 10~ 8 . The numerical results are listed in 
Tables / to Table VI for different values of n, m and different eccentricities of some members of the 
solar system. In these tables 8 2 A {5 2 B ) represents the sum of the squares of the residuals of Equation (19) 
for A's(B's) coefficients, a coe f f.A{vcoef f.s) represents the common standard error of Equation( 21) 
for A's(B's) coefficients, finally Qa{Qb) represents the average squared distance between the exact 
and least- squares values of Equation (23) for A's(B's) coefficients. 

4 CONCLUSION 

In concluding the present paper, pure numerical method is developed for computing Hansen coeffi- 
cients by using recursive harmonic analysis technique. The precision criteria which are: the variance of 
the fit, the standard errors of the coefficients and the average squared distance between the exact and 
least squares values, are all very satisfactory. The method is not only provide materials for computing 
Hansen's and also Hansen's like expansions but also can be used to check the accuracy of the different 
algorithms that already exist. 
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Table I: Hansen Coefficients for The Planet Earth : 
e = .016708617, n = -3, m = 6 



k 


A k 


B k 





-2.80505 x 10~ iu 




1 


-1.69927 x 10 -10 


-1.69927 x 10" 


2 


1.31491 x 10 -7 


1.31491 x 10~ 7 


3 


-0.0000259261 


-0.0000259261 


4 


0.00209013 


0.00209013 


5 


-0.0749101 


-0.0749101 


6 


0.99039 


0.99039 


7 


0.124591 


0.124591 


8 


0.00917108 


0.00917108 


9 


0.000516607 


0.000516607 


10 


0.0000246565 


0.0000246565 


1 1 


1.05004 x 10" 6 


1.05004 x 10" 6 



d\ = 8.52652 x 10~ 14 d B = 6.39488 x 10~ 14 

Soo ee .A = 4.37729 x 10~ 9 s xee . B = 3.79085 x 10~ 9 

Qa = 2.10768 x 10~ 16 Q B = 1.8076 x 10~ 16 



Table II: Hansen Coefficients for The Planet Pluto : 
e = 0.249050, n = -3, m = 6 



k 


A k 


B k 





0.0508079 




1 


-0.325005 


-0.319177 


2 


0.969155 


0.969203 


3 


-1.34716 


-1.34716 


4 


0.536896 


0.536896 


5 


0.248699 


0.248699 


6 


0.0765758 


0.0765758 


7 


0.0211903 


0.0211903 


8 


0.0056458 


0.0056458 


9 


0.00148234 


0.00148234 


10 


0.000386931 


0.000386931 


11 


0.000100722 


0.000100722 


12 


0.0000261571 


0.0000261571 


13 


6.76933 x 10" 6 


6.76933 x 10"' 



d\ = 1.62174 x 10~ iu d B = 1.6226 x 10~ ilJ 

Sooee.A = 1.93084 X 10" 7 Sooee.S = 1.93135 X 10" 7 

Qa = 4.84659 x 10" 13 Q B = 4.84914 x 10~ 13 
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Table III: Hansen Coefficients for The Asteroid Ceres: 
e = 0.078, n = 8, m = 2 



B k 



0.0854431 

-0.492936 

1.08609 

-0.157994 

0.00140598 

0.000192711 

0.0000113508 



-0.479094 

1.08564 

-0.157993 

0.00140603 

0.000192714 

0.0000113509 



d\ = 4.26326 x 10~ 14 
Sooee.A = 3.02792 x 10~ £ 
Qa = 6.41781 x 10~ 17 



d B = 4.26326 x 10~ i4 
Sooee.S = 3.02792 x 10" 
Qb = 6.41781 x 10~ 17 



Table IV: Hansen Coefficients for The Asteroid Sekhmet : 
e = 0.296, n = -l,m = 5 



k 


A h 


Bk 





-0.0000795273 




1 


0.00983893 


0.00983416 


2 


-0.114213 


-0.114214 


3 


0.431088 


0.431088 


4 


-0.482649 


-0.482649 


5 


-0.260258 


-0.260258 


6 


0.191795 


0.191795 


7 


0.410314 


0.410314 


8 


0.411965 


0.411965 


9 


0.318733 


0.318733 


10 


0.213837 


0.213837 


11 


0.130854 


0.130854 


12 


0.0750457 


0.0750457 


13 


0.0410069 


0.0410069 


14 


0.021582 


0.021582 


15 


0.0110231 


0.0110231 


16 


0.00549385 


0.00549385 


17 


0.00268279 


0.00268279 


18 


0.00128767 


0.00128767 


19 


0.000608982 


0.000608982 


20 


0.000284349 


0.000284349 


21 


0.000131294 


0.000131294 


22 


0.0000600289 


0.0000600289 


23 


0.0000272069 


0.0000272069 


24 


0.0000122351 


0.0000122351 


25 


5.46368 x 10~ 6 


5.46368 x 10" 



d\ = 3.64729 x 10" ilJ d B = 3.64665 x 10" lu 

Sooee.A = 3.11867 X 10~ 7 Socee.B = 3.1184 X HP 7 

Qa = 2.43152 x 10~ 12 Q B = 2.4311 x 10~ 12 
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Table V: Hansen Coefficients for The Comet Wild2 : 
e = 0.541, n = 3,m = 2 



fc 


A k 


B k 





-0.187235 




1 


0.954443 


0.943797 


2 


-1.75935 


-1.75982 


3 


1.04451 


1.04448 


4 


0.271626 


0.271625 


5 


-0.101094 


-0.101092 


6 


-0.158792 


-0.158791 


7 


-0.112669 


-0.112669 


8 


-0.05443 


-0.0544298 


9 


-0.0112087 


-0.0112086 


10 


0.0142725 


0.0142726 


1 1 


0.026263 


0.026263 


12 


0.0296931 


0.0296931 


13 


0.0283862 


0.0283862 


14 


0.0248699 


0.0248699 


15 


0.0206484 


0.0206484 


16 


0.0165297 


0.0165297 


17 


0.0128893 


0.0128893 


18 


0.00985376 


0.00985376 


19 


0.00741822 


0.00741822 


20 


0.005551671 


0.005551671 


21 


0.00406198 


0.00406198 


22 


0.00296636 


0.00296636 


23 


0.00215139 


0.00215139 


24 


0.00155123 


0.00155123 


25 


0.00111291 


0.00111291 


26 


0.000795003 


0.000795003 


27 


0.000565771 


0.000565771 


28 


0.000401309 


0.000401309 


29 


0.000283824 


0.000283824 


30 


0.000200214 


0.000200214 


31 


0.000140908 


0.000140908 


32 


0.0000989644 


0.0000989644 


33 


0.0000693761 


0.0000693761 


34 


0.000048552 


0.000048552 


35 


0.0000339265 


0.0000339265 


36 


0.0000236736 


0.0000236736 


37 


0.0000164982 


0.0000164982 


38 


0.0000114842 


0.0000114842 


39 


7.9854 x 10 -6 


7.9854 x 10~ 6 



d\ = 1.18559 x 10~ s d\ = 1.18562 x 10~ s 

Sooee.A = 4.05228 X 10" 7 Socee.B = 4.05232 x 10~ 7 

Qa = 6.40418 x IP" 12 Q B = 6.4043 x 10" 12 
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Table VI: Hansen Coefficients for The Comet Lexell : 
e = 0.786, n = 8, m = 4 



k 


A k 


B k 





28.4068 




1 


-47.0631 


-25.693 


2 


23.9162 


21.1464 


3 


-4.70405 


-4.84203 


4 


-0.605464 


-0.619241 


5 


-0.0262285 


-0.0283248 


6 


0.0293643 


0.0289445 


7 


0.0217703 


0.0216686 


8 


0.0121887 


0.0121605 


9 


0.00637396 


0.00636531 


10 


0.00325957 


0.00325672 


11 


0.00164622 


0.00164524 


12 


0.000817397 


0.000817052 


13 


0.000392683 


0.000392564 


14 


0.000176143 


0.000176105 


15 


0.0000672271 


0.0000672181 


16 


0.0000140575 


0.0000140581 


17 


-0.0000103329 


-0.0000103297 


18 


-0.0000200483 


-0.000020045 


19 


-0.0000224845 


-0.0000224818 


20 


-0.0000215081 


-0.0000215061 


21 


-0.0000191168 


-0.0000191154 


22 


-0.0000163167 


-0.0000163158 


23 


-0.0000135893 


-0.0000135887 


24 


-0.0000111418 


-0.0000111414 


25 


-9.0409 x 10~ 6 


-9.04071 x 10" 



d A = 7.42148 x 10" a d z B = 7.33417 x 10" a 

Sooee.A = 1.40679 x 10" 6 Sacee . B = 1.39849 x 10" 

Qa = 4.94765 x IP' 11 Qb — 4.88944 x IP" 11 



